{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Example: Virtual circuit calculations\n", "\n", "---\n", "\n", "This example notebook will demonstrate how to construct **virtual circuits** (VCs) for a **MAST-U-like tokamak** equilibrium.\n", "\n", "### What are Virtual Circuits and why are they useful?\n", "\n", "VCs can be used to identify which (active) poloidal field (PF) coils have the most significant impact on a set of specified plasma shape parameters (which we will refer to henceforth as **targets**). The targets are equilibrium-related quantities such as the inner or outer midplane radii $R_{in}$ or $R_{out}$ (more will be listed later on). \n", "\n", "More formally, the **virtual circuit** (VC) matrix $V$ for a given equilibrium (and chosen set of shape targets and PF coils) is defined as\n", "\n", "$$ V = (S^T S)^{-1} S^T, $$\n", "\n", "where $S$ is the **shape** (Jacobian) matrix:\n", "\n", "$$ S_{i,j} = \\frac{\\partial T_i}{\\partial I_j}. $$\n", "\n", "Here, $T_i$ is the $i$ th target and $I_j$ is the current in the $j$ th PF coil. We note that $V$ is simply the Moore-Penrose pseudo-inverse of $S$. In FreeGSNKE, $S$ will be calculated using finite difference (see below). \n", "\n", "### How can these VCs be used?\n", "\n", "Once we know the VC matrix $V$ for a given equilibrium (and its associated target values $\\vec{T}$), we can specify a perturbation in the targets $\\vec{\\Delta T}$ and calculate the change in coil currents required to acheive the new targets. The shifts in coil currents can be found via:\n", "\n", "$$ \\vec{\\Delta I} = V \\vec{\\Delta T}. $$\n", "\n", "Using $\\vec{\\Delta I}$ we can perturb the coil currents in our equilibrium, re-solve the static forward Grad-Shafranov (GS) problem, and observe how the targets (call them $\\vec{T}_{new}$) have changed \n", "in the new equilibrium vs. the old targets $\\vec{T}$.\n", "\n", " ---\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Generate a starting equilibrums\n", "\n", "Firstly, we need an equilbirium to test the VCs on." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "# build machine\n", "from freegsnke import build_machine\n", "tokamak = build_machine.tokamak(\n", " active_coils_path=f\"../machine_configs/MAST-U/MAST-U_like_active_coils.pickle\",\n", " passive_coils_path=f\"../machine_configs/MAST-U/MAST-U_like_passive_coils.pickle\",\n", " limiter_path=f\"../machine_configs/MAST-U/MAST-U_like_limiter.pickle\",\n", " wall_path=f\"../machine_configs/MAST-U/MAST-U_like_wall.pickle\",\n", " magnetic_probe_path=f\"../machine_configs/MAST-U/MAST-U_like_magnetic_probes.pickle\",\n", ")\n", "\n", "# initialise the equilibrium\n", "from freegsnke import equilibrium_update\n", "eq = equilibrium_update.Equilibrium(\n", " tokamak=tokamak,\n", " Rmin=0.1, Rmax=2.0, # Radial range\n", " Zmin=-2.2, Zmax=2.2, # Vertical range\n", " nx=65, # Number of grid points in the radial direction (needs to be of the form (2**n + 1) with n being an integer)\n", " ny=129, # Number of grid points in the vertical direction (needs to be of the form (2**n + 1) with n being an integer)\n", " # psi=plasma_psi\n", ") \n", "\n", "# initialise the profiles\n", "from freegsnke.jtor_update import ConstrainPaxisIp\n", "profiles = ConstrainPaxisIp(\n", " eq=eq, # equilibrium object\n", " paxis=8e3, # profile object\n", " Ip=6e5, # plasma current\n", " fvac=0.5, # fvac = rB_{tor}\n", " alpha_m=1.8, # profile function parameter\n", " alpha_n=1.2 # profile function parameter\n", ")\n", "\n", "# load the nonlinear solver\n", "from freegsnke import GSstaticsolver\n", "GSStaticSolver = GSstaticsolver.NKGSsolver(eq) \n", "\n", "# set the coil currents\n", "import pickle\n", "with open('data/simple_diverted_currents_PaxisIp.pk', 'rb') as f:\n", " current_values = pickle.load(f)\n", "for key in current_values.keys():\n", " eq.tokamak.set_coil_current(key, current_values[key])\n", "\n", "# carry out the foward solve to find the equilibrium\n", "GSStaticSolver.solve(eq=eq, \n", " profiles=profiles, \n", " constrain=None, \n", " target_relative_tolerance=1e-4)\n", "\n", "# updates the plasma_psi (for later on)\n", "eq._updatePlasmaPsi(eq.plasma_psi)\n", "\n", "# plot the resulting equilbria \n", "fig1, ax1 = plt.subplots(1, 1, figsize=(4, 8), dpi=80)\n", "ax1.grid(True, which='both')\n", "eq.plot(axis=ax1, show=False)\n", "eq.tokamak.plot(axis=ax1, show=False)\n", "ax1.set_xlim(0.1, 2.15)\n", "ax1.set_ylim(-2.25, 2.25)\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Initialise the VC class\n", "\n", "We initialise the VC class and tell it to use the static solver we already initialised above. This will be used to repeatedly (and rapidly) to solve the static GS problem when calculating the finite differences for the shape matrix." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from freegsnke import virtual_circuits\n", "\n", "VCs = virtual_circuits.VirtualCircuitHandling()\n", "VCs.define_solver(GSStaticSolver)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Define shape targets\n", "\n", "Next we need to define the shape targets (i.e. our quantities of interest) that we wish to monitor. \n", "\n", "The way we do this is exactly the same as was seen in Example05b with the \"plasma_descriptors\" function. The user just needs to define a function that takes in the `eq` as a (single) input and returns an array of the shape targets (sometimes called parameters) of interest. See the function defined in the next cell." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# define the descriptors function (it should return an array of values and take in an eq object)\n", "def plasma_descriptors(eq):\n", "\n", " # inboard/outboard midplane radii\n", " RinRout = eq.innerOuterSeparatrix()\n", "\n", " # inboard midplane radius\n", " Rout = eq.innerOuterSeparatrix()[0]\n", "\n", " # find lower X-point\n", " # define a \"box\" in which to search for the lower X-point\n", " XPT_BOX = [[0.33, -0.88], [0.95, -1.38]]\n", "\n", " # mask those points\n", " xpt_mask = (\n", " (eq.xpt[:, 0] >= XPT_BOX[0][0])\n", " & (eq.xpt[:, 0] <= XPT_BOX[1][0])\n", " & (eq.xpt[:, 1] <= XPT_BOX[0][1])\n", " & (eq.xpt[:, 1] >= XPT_BOX[1][1])\n", " )\n", " xpts = eq.xpt[xpt_mask, 0:2].squeeze()\n", " if xpts.ndim > 1 and xpts.shape[0] > 1:\n", " opt = eq.opt[0, 0:2]\n", " dists = np.linalg.norm(xpts - opt, axis=1)\n", " idx = np.argmin(dists) # index of closest point\n", " Rx, Zx = xpts[idx, :]\n", " else:\n", " Rx, Zx = xpts\n", "\n", " return np.array([RinRout[0], RinRout[1], Rx, Zx])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Assign these shape targets some decriptive names (and ensure the names are in the correct order!). You can then call the function to test it works correctly. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "target_names = [\"Rin\", \"Rout\", \"Rx\", \"Zx\"]\n", "target_values = plasma_descriptors(eq)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us visualise the targets on our equilibrium. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# plot the resulting equilbria \n", "fig1, ax1 = plt.subplots(1, 1, figsize=(4, 8), dpi=80)\n", "ax1.grid(True, which='both')\n", "ax1.contour(eq.R, eq.Z, eq.psi(), levels=[eq.psi_bndry], colors='r')\n", "eq.tokamak.plot(axis=ax1, show=False)\n", "ax1.plot(eq.tokamak.wall.R, eq.tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", "\n", "ax1.scatter(target_values[0], 0.0, s=100, color='green', marker='x', zorder=20, label=target_names[0])\n", "ax1.scatter(target_values[1], 0.0, s=100, color='blue', marker='x', zorder=20, label=target_names[1])\n", "ax1.scatter(target_values[2], target_values[3], s=100, color='m', marker='*', zorder=20, label=f\"({target_names[2]},{target_names[3]})\")\n", "\n", "ax1.set_xlim(0.1, 2.15)\n", "ax1.set_ylim(-2.25, 2.25)\n", "ax1.legend(loc=\"upper right\")\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Calculating the VCs\n", "\n", "Now we've defined the targets we're interested in, we can begin calculating the **shape** and **virtual circuit** matrices. The following is a brief outline of how they're calculated:\n", "\n", "##### 1. Initialise and solve the base equilibrium\n", "Before computing derivatives, the static forward GS solver is run with some initial coil currents. Following this we store the:\n", "- equilibrium state `eq` and plasma profiles `profiles`.\n", "- plasma current vector $I_y$ (which defines the amount of plasma current at each computational grid point, restricted to the limiter region to save computational resources).\n", "- target values $\\vec{T} = [T_1,\\ldots,T_{N_T}]$ are evaluated.\n", "\n", "This establishes a reference state before we perturb the coil currents - we have already carried out these steps above (the $I_y$ vector is stored within `profiles`).\n", "\n", "##### 2. Find the appropriate coil current perturbations\n", "Next, we need to identify the appropriate coil current perturbations $\\vec{\\delta I}$ to use in the finite differences. To do this, we aim to scale the starting guess `starting_dI` (denoted $\\vec{\\delta I}^{start}$) according to the target tolerance on the plasma current vector (`target_dIy` = $\\delta I_y^{target}$). While the user can choose `starting_dI`, the default is given by\n", "\n", "$$ \\vec{\\delta I}^{start} := | \\vec{I} | \\times \\delta I_y^{target}. $$\n", "\n", "For each coil current $j \\in \\{1,\\ldots,N_c \\}$, this starting guess is scaled as follows:\n", "1. Perturb coil current $j$: \n", "\n", "$$ I_j^{new} := I_j + \\delta I_j^{start}.$$\n", "\n", "2. Solve the equilibrium with the updated $j$ th coil current (others are left unchanged) and store the plasma current vector $\\vec{I}_y^{new}$.\n", "\n", "3. The starting guess for the $j$ th coil current perturbation is then scaled by the relative norm of the plasma current change (and the predefined target tolerance `target_dIy`):\n", "\n", "$$ \\delta I_j^{final} = \\delta I_j^{start} \\times \\left( \\delta I_y^{target} \\frac{\\| \\vec{I}_y^{start} \\|}{\\| \\vec{I}_y^{new} - \\vec{I}_y^{start} \\|} \\right).$$\n", "\n", "If this relative norm is larger than $\\delta I_y^{target}$, then the starting guess $\\delta I_j^{start}$ needs to be made smaller (and vice versa). \n", "\n", "After this, we have our scaled perturbations $\\vec{\\delta I}^{final}$ ready to use in the finite differences. \n", "\n", "##### 3. Find the finite differences (i.e. the shape matrix)\n", "For each coil current $j \\in \\{1,\\ldots,N_c \\}$:\n", "\n", "1. Perturb coil current $j$: \n", "\n", "$$ I_j^{new} := I_j + \\delta I_j^{final}.$$\n", "\n", "2. Solve the equilibrium with the updated $j$ th coil current (others are left unchanged) and store the plasma current vector $\\vec{I}_y^{new}$.\n", "\n", "3. Using the new target values $\\vec{T}^{new}$ from the equilibrium, calculate the $j$ th column of the shape matrix:\n", "\n", "$$ S_{:,j} = \\frac{\\vec{T}^{new} - \\vec{T}}{\\delta I_j^{final}}. $$\n", "\n", "This gives the sensitivity of each of the targets $T_i$ with respect to a change in the $j$ th coil current.\n", "\n", "Note that we can also obtain the Jacobian matrix of the plasma current vector (using $\\vec{I}_y^{new}$) with respect to the coil currents: $\\frac{\\partial \\vec{I}_y}{\\partial \\vec{I}}$.\n", "\n", "##### 4. Find the virtual circuit matrix\n", "Once the full shape matrix $S \\in \\Reals^{N_T \\times N_c} $ is known, the **virtual circuit matrix** is computed as:\n", "\n", "$$ V = (S^T S)^{-1} S^T \\in \\Reals^{N_c \\times N_T}.$$\n", "\n", "This matrix provides a mapping from requested shifts in the targets to the shifts in the coil currents required." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# next we need to define which coils we wish to calculate the shape (finite difference) derivatives (and therefore VCs) for\n", "print(\"All active coils :\", eq.tokamak.coils_list[0:12])\n", "\n", "# here we'll look at a subset of the coils and the targets\n", "coils = ['PX', 'D1', 'D2', 'D3', 'Dp', 'D5', 'D6', 'D7', 'P4', 'P5']\n", "print(\"Coils to calculate VCs for :\", coils)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# calculate the shape and VC matrices using the in-built method\n", "VCs.calculate_VC(\n", " eq=eq,\n", " profiles=profiles,\n", " coils=coils,\n", " target_names=target_names,\n", " target_calculator=plasma_descriptors,\n", " starting_dI=None, \n", " min_starting_dI=50,\n", " verbose=True,\n", " name=\"VC_for_lower_targets\", # name for the VC\n", " )" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# these are the finite difference derivatives of the targets wrt the coil currents\n", "shape_matrix = 1.0*VCs.VC_for_lower_targets.shape_matrix\n", "print(f\"Dimension of the shape matrix is {shape_matrix.shape} [no. of targets x no. of coils]. Has units [m/A]\")\n", "\n", "# these are the VCs corresponding to the above shape matrix\n", "VCs_matrix = 1.0*VCs.VC_for_lower_targets.VCs_matrix\n", "print(f\"Dimension of the virtual circuit matrix is {VCs_matrix.shape} [no. of coils x no. of targets]. Has units [A/m].\")\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### How do we make use of the VCs?\n", "\n", "Now that we have the VCs, we can use them to idenitfy the coil current shifts required to change the targets by a certain amount. \n", "\n", "For example, we will ask for shifts in a few shape targets and observe what happens to the equilibrium once we apply the new currents from the VCs. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# let's remind ourselves of the targets\n", "print(target_names)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# let's set the requested shifts (units are metres) for the full set of targets\n", "requested_target_shifts = [0.05, -0.05, 0.0, 0.0]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We specify a perturbation in the targets $\\vec{\\Delta T}$ (the shifts above) and calculate the change in coil currents required to achieve this by calculating:\n", "\n", "$$ \\vec{\\Delta I} = V \\vec{\\Delta T}. $$\n", "\n", "Using $\\vec{\\Delta I}$ we then perturb the coil currents in our original equilibrium, re-solve the static forward Grad-Shafranov (GS) problem, and return. This is all done in the following cell. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# we apply the VCs to the desired shifts, using the VC_object we just created\n", "eq_new, profiles_new, new_target_values, old_target_values = VCs.apply_VC(\n", " eq=eq,\n", " profiles=profiles,\n", " VC_object=VCs.VC_for_lower_targets,\n", " requested_target_shifts=requested_target_shifts,\n", " verbose=True,\n", ")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# plot the resulting equilbria \n", "fig1, ax1 = plt.subplots(1, 1, figsize=(4, 8), dpi=80)\n", "ax1.grid(True, which='both')\n", "ax1.plot(eq.tokamak.wall.R, eq.tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", "eq.tokamak.plot(axis=ax1, show=False)\n", "\n", "ax1.contour(eq.R, eq.Z, eq.psi(), levels=[eq.psi_bndry], colors='r')\n", "ax1.contour(eq.R, eq.Z, eq_new.psi(), levels=[eq_new.psi_bndry], colors='b')\n", "\n", "ax1.set_xlim(0.1, 2.15)\n", "ax1.set_ylim(-2.25, 2.25)\n", "\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can then measure the accuracy of the VCs by observing the difference between the requested change in shape targets vs. those enacted by the VCs. \n", " \n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "for i in range(len(target_names)):\n", " print(f\"Difference in {target_names[i]} = {np.round(new_target_values[i] - old_target_values[i],3)} vs. requested = {(requested_target_shifts)[i]}.\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the plots below, we can see the actual shifts by the VCs are almost exactly as those requested (for the targets with actual shifts). Note how the unshifted targets also shift under the VCs due to the nonlinear coupling between targets and coil currents. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# plots\n", "rel_diff = np.abs(np.array(new_target_values) - np.array(old_target_values))/np.abs(old_target_values)\n", "\n", "fig1, ((ax1), (ax2)) = plt.subplots(2, 1, figsize=(12,12), dpi=80)\n", "\n", "ax1.grid(True, which='both', alpha=0.5)\n", "ax1.scatter(target_names, requested_target_shifts, color='red', marker='o', s=150, label=\"Requested\")\n", "ax1.scatter(target_names, np.array(new_target_values) - np.array(old_target_values), color='royalblue', marker='o', s=75, label=\"Actual\")\n", "ax1.set_xlabel(\"Target\")\n", "ax1.set_ylabel(\"Shift [m]\")\n", "ax1.legend()\n", "# ax1.set_ylim([-max(targets_shift+non_standard_targets_shift)*1.1, max(targets_shift+non_standard_targets_shift)*1.1])\n", "\n", "ax2.grid(True, which='both', alpha=0.5)\n", "ax2.scatter(target_names, rel_diff, color='red', marker='o', s=150, edgecolors='black', label=\"Requested\")\n", "ax2.set_xlabel(\"Target\")\n", "ax2.set_ylabel(\"Relative shift\")\n", "labels = ax2.get_xticklabels()\n", "ax2.set_yscale('log')\n", "ax2.set_ylim([1e-6, 1e-0])\n", "\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Alternative VCs\n", "\n", "Note that when the `VC_calculate` method is called, a `name` can be provided. This will store the key data in a subclass with `name`. This useful for recalling from which targets and coils a given shape matrix or VC matrix was calculated. \n", "\n", "Note that if no `name` is provided the default is used (\"latest_VC\")." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# display the attributes\n", "print(VCs.VC_for_lower_targets.__dict__.keys())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can then follow up by calculating alternative VCs for different target and coil combinations. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# define the descriptors function (it should return an array of values and take in an eq object)\n", "def plasma_descriptors_new(eq):\n", "\n", " # inboard/outboard midplane radii\n", " RinRout = eq.innerOuterSeparatrix()\n", "\n", " # find lower X-point\n", " # define a \"box\" in which to search for the lower X-point\n", " XPT_BOX = [[0.33, -0.88], [0.95, -1.38]]\n", "\n", " # mask those points\n", " xpt_mask = (\n", " (eq.xpt[:, 0] >= XPT_BOX[0][0])\n", " & (eq.xpt[:, 0] <= XPT_BOX[1][0])\n", " & (eq.xpt[:, 1] <= XPT_BOX[0][1])\n", " & (eq.xpt[:, 1] >= XPT_BOX[1][1])\n", " )\n", " xpts = eq.xpt[xpt_mask, 0:2].squeeze()\n", " if xpts.ndim > 1 and xpts.shape[0] > 1:\n", " opt = eq.opt[0, 0:2]\n", " dists = np.linalg.norm(xpts - opt, axis=1)\n", " idx = np.argmin(dists) # index of closest point\n", " Rx, Zx = xpts[idx, :]\n", " else:\n", " Rx, Zx = xpts\n", "\n", " return np.array([RinRout[0], Rx, Zx])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# this time we use only the divertor coils and a subset of targets\n", "coils = [\"PX\", 'D1', 'D2', 'D3', 'Dp', 'D5', 'D6', 'D7']\n", "target_names =['Rin', 'Rx', 'Zx']\n", "\n", "\n", "# calculate the shape and VC matrices as follows\n", "VCs.calculate_VC(\n", " eq=eq,\n", " profiles=profiles,\n", " coils=coils,\n", " target_names=target_names,\n", " target_calculator=plasma_descriptors_new,\n", " starting_dI=None,\n", " min_starting_dI=50,\n", " verbose=True,\n", " name=\"VC_for_upper_targets\", # name for the VC\n", " )" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# display the attributes\n", "print(VCs.VC_for_upper_targets.__dict__.keys())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can now apply this VC using the `apply_VC` method." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# we'll shift the upper strikepoint (R value) \n", "requested_target_shifts = [0.0, -0.025, 0.03]\n", "\n", "# we apply the VCs to the some desired shifts, using the VC_object we just created\n", "eq_new, profiles_new, new_target_values, old_target_values = VCs.apply_VC(\n", " eq=eq,\n", " profiles=profiles,\n", " VC_object=VCs.VC_for_upper_targets,\n", " requested_target_shifts=requested_target_shifts,\n", " verbose=True,\n", ")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# plot the resulting equilbria \n", "fig1, ax1 = plt.subplots(1, 1, figsize=(4, 8), dpi=80)\n", "ax1.grid(True, which='both')\n", "ax1.plot(eq.tokamak.wall.R, eq.tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", "eq.tokamak.plot(axis=ax1, show=False)\n", "\n", "ax1.contour(eq.R, eq.Z, eq.psi(), levels=[eq.psi_bndry], colors='r')\n", "ax1.contour(eq.R, eq.Z, eq_new.psi(), levels=[eq_new.psi_bndry], colors='b')\n", "\n", "ax1.set_xlim(0.1, 2.15)\n", "ax1.set_ylim(-2.25, 2.25)\n", "\n", "plt.tight_layout()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# plots\n", "rel_diff = np.abs(np.array(new_target_values) - np.array(old_target_values))/np.abs(old_target_values)\n", "\n", "fig1, ((ax1), (ax2)) = plt.subplots(2, 1, figsize=(12,12), dpi=80)\n", "\n", "ax1.grid(True, which='both', alpha=0.5)\n", "ax1.scatter(target_names, requested_target_shifts, color='red', marker='o', s=150, label=\"Requested\")\n", "ax1.scatter(target_names, np.array(new_target_values) - np.array(old_target_values), color='royalblue', marker='o', s=75, label=\"Actual\")\n", "ax1.set_xlabel(\"Target\")\n", "ax1.set_ylabel(\"Shift [m]\")\n", "ax1.legend()\n", "# ax1.set_ylim([-max(targets_shift+non_standard_targets_shift)*1.1, max(targets_shift+non_standard_targets_shift)*1.1])\n", "\n", "ax2.grid(True, which='both', alpha=0.5)\n", "ax2.scatter(target_names, rel_diff, color='red', marker='o', s=150, edgecolors='black', label=\"Requested\")\n", "ax2.set_xlabel(\"Target\")\n", "ax2.set_ylabel(\"Relative shift\")\n", "labels = ax2.get_xticklabels()\n", "ax2.set_yscale('log')\n", "ax2.set_ylim([1e-6, 1e-0])\n", "\n", "plt.tight_layout()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "freegsnke4e", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.16" } }, "nbformat": 4, "nbformat_minor": 4 }